Spinwave damping in the two-dimensional ferromagnetic XY model 
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The effect of damping of spinwaves in a two-dimensional 
classical ferromagnetic XY model is considered. The damping 
rate T q is calculated using the leading diagrams due to the 
quartic-order deviations from the harmonic spin Hamiltonian. 
The resulting four-dimensional integrals are evaluated by ex- 
tending the techniques developed by Gilat and others for spec- 
tral density types of integrals. F q is included into the memory 
function formalism due to Reiter and Solander, and Menezes, 
to determine the dynamic structure function S(q, u). For the 
infinite sized system, the memory function approach is found 
to give non-divergent spinwave peaks, and a smooth nonzero 
background intensity ("plateau" or distributed intensity) for 
the whole range of frequencies below the spinwave peak. The 
background amplitude relative to the spinwave peak rises with 
temperature, and eventually becomes higher than the spin- 
wave peak, where it appears as a central peak. For finite- 
sized systems, there are multiple sequences of weak peaks on 
both sides of the spinwave peaks whose number and positions 
depend on the system size and wavevector in integer units 
of 27r/L . These dynamical finite size effects are explained in 
the memory function analysis as due to either spinwave dif- 
ference processes below the spinwave peak or sum processes 
above the spinwave peak. These features are also found in 
classical Monte Carlo - Spin-Dynamics simulations. 

PACS numbers: 75.10.Hk,75.30.Ds,75.40.Gb,75.40.Mg 



I. INTRODUCTION 

The two-dimensional (2D) ferromagnetic XY model 
has had a long history of studies concerning the dynamics 
of spinwafies at low temperatures. It was first analyzed 
by VillainEJ and by Moussa and Villainl, for low temper- 
atures, in the harmonic approximation. They found that 
the dynamic structure function S(q, u>) has a spinwave 
peak of the form 



1-1+V/2 



(1) 



where r\ is the critical exponent describing the decay of 
the static spin-spin correlation function, and cj q is t 
spinwave frequency at wavevector q. Nelson and Fishc 
treated the model in a fixed-length hydrodynamic de- 
scription, obtaining the following expression around the 
spinwave peak, 
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where for small q, ujq = cq, with c being the spinwave 
velocity. These theoretical approaches led to spinwave 
peaks diverging unphysically at the spinwave frequency. 
Here we describe an alternative theoretical description in 
which the divergence does not occur. 

At higher temperatures, the spinwave or harmonic ap- 
proximation begins to fail, and XY magnets are known to 
undergo a topological phase transitiorofl associated with 
the unbinding of vortices. We generally consider temper- 
atures below this unbinding temperature Tbkt, where a 
spinwave approximation will be valid. At very low tem- 
peratures any nonlinear effects should be very weak, do 
not lead to vortex production, and can be included as 
the damping or scattering of the spinwave modes, that 
results in a spinwave lifetime. Thus our approach will 
be to start from the harmonic approximation to the spin 
Hamiltonian, Hq, and then include the damping that re- 
sults from the terms fourth order in spin operators (Hi) 
in an expansion of the original full spin Hamiltonian, H. 
The resulting damping strength T q will imply a finite 
lifetime r = T" 1 for the spinwave modes. We are inter- 
ested to investigate whether the presence of such lifetime 
effects can eliminate the divergence of the dynamic struc- 
ture function S (q, u>) at the spinwave frequency. 

We calculate the spinwave damping by a Feynmann 
diagram technique for finite temperatures, applied to the 
fourth-order Hamiltonian Hi , using the leading diagrams 
to second order in H\. The only diagrams we include 
are those that give a contribution to the imaginary part 
of the spinwave energy, i.e., to the damping or inverse 
lifetime. These are the terms that involve a dynami- 
cal effect on the spinwave energies. The associated inte- 
gral expressions for these diagrams require the evaluation 
of four-dimensional integrals involving a delta-function. 
Methods for evaluating this type of integral in three di- 
mensions have been developed by Gilat and co-workers; 
here we show how to extend their work and apply the 
same technique in four dimensions. 

In principle, we could also use the Feynmann tech- 
nique to determine the leading changes in the real part 
of the spinwave energy, Aw q , due to the perturbation 
H\. This real part would result as the usual effect due 
to finite temperature: a shift in the spinwave frequen- 
cies to lower values as the temperature increases. It 
is essentially a static effect, primarily due to the renor- 
malization of the nearest neighbor exchange interaction 
strength J with temperature. Calculation of these spin- 
wave frequency shifts by the Feynmann technique would 
require including diagrams first order in H\ ; the so-called 
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single-loop diagrams. In this paper, however, we are pri- 
marily concerned with accounting for dynamical effects. 
We will include the renormalization of the spinwave en- 
ergies using the memory function formalism, combined 
with some input information from classical Monte Carlo 
(MC) data. In the memory function approach, the posi- 
tion of the spinwave peak in S"(q, u)) does not occur ex- 
actly at the zero-temperature dispersion relation, w q , but 
rather, close to the value J (w 2 ) q < w q , where (w 2 ) q is 
the temperature-dependent second moment. (w 2 ) q can 
be estimated from static quantities, either using low- 
temperature analytic expressions or results from Monte 
Carlo, the latter being necessary for temperatures ap- 
proaching Tbkt- The memory function approach already 
includes to leading order the exchange strength renormal- 
ization, as can be seen in the shift in the position of the 
spinwave peaks. Thus we use it rather than making any 
renormalization of the exchange J. 

A few comments can be made about how the damp- 
ing T q is included into the memory function calculation. 
Note that the damping calculation sketched above is com- 
pletely independent of the memory function calculation; 
the latter is a formalism for calculating dynamic correla- 
tions for a given Hamiltonian and operators. Due to the 
XY symmetry of the model, the natural coordinates of 
the spins are the in-plane angle <fi and out-of-plane spin 
component S z . Fourier transforming and using a stan- 
dard canonical transformation, the Hamiltonian is writ- 
ten in terms of creation and annihilation operators and 
a q , with time dependencies varying as e IWqt and e _!IJqi , 
respectively. Once damping is included, the time depen- 
dencies become e *( w q+ ir q)* arLC j e — i(w q — ir,,)^ reS p ec ti V ely. 
Henceforth the calculation proceeds as usual in the mem- 
ory function approach, except for this small modification 
of the operator time dependencies via the shifted opera- 
tor frequencies, w q — > w q ± iF q . 

In the memory function calculation, spinwave sum and 
difference processes enter. These are characterized by 
combinations of frequencies, fi± = wa +k ± cjq_ k , where 
q is the wavevector at which the dynamic correlations are 
measured, and k is summed over to include all possible 
such processes. Based on the above modifications of the 
operator frequencies, there also results the modification 
for these processes, £2± — * Q± + i(ra +k + Tq_ k ). The 
damping rates for different wavevectors always add for 
both types of processes. 

The paper is organized as follows: First, in Sec. |J, 
we present the model and develop the expansion of the 
Hamiltonian to quartic order in creation/annihilation op- 
erators. In Sec. E|, the expressions for the damping via 
the Feynmann technique are presented. Then we show 
how these are evaluated using the Gilat procedure. In 
Sec. IV, the key aspects of the memory function calcu- 



lation of S(q, u>) are reviewed, together with a focus on 
the modification due to the damping. Some details of the 
Monte Carlo and spin dynamics simu latio ns are summa- 
rized in Sec. [v|. This is followed in Sec. VI by presentation 



of the results both for very large and relatively smaller 
systems, in order to discuss the dynamical finite size ef- 
fects. The memory function results also are compared 
with standard Monte Carlo - spin-dynamics calculations 
of the dynamic correlation functions in order to evaluate 
the validity and accu racy of the method. Our conclusions 



are presented in Sec. VII 



II. THE MODEL 

The ferromagnetic spin Hamiltonian under considera- 
tion is 



where the first sum is over all sites n of a 2D square 
lattice, the second is over all nearest neighbor displace- 
ments a, J is the nearest neighbor coupling strength, 
and each S n is a three-component classical spin-vector. 
Because each bond is counted twice, a factor of \ is in- 
cluded on J. We primarily concentrate on the XY model, 
with anisotropy parameter A = 0, however, for complete- 
ness, formulas covering the easy-plane anisotropy range, 
< A < 1, are presented. 



A. Expansion of Hamiltonian 

In order to develop the damping as a perturbation of 
the harmonic approximation, we need to expand the full 
spin Hamiltonian, including the harmonic (quadratic) 
and quartic deviations from the ground state. The 
ground state corresponds to spins aligned in any direction 
within the XY plane. Writing the spins in Hamiltonian 
(0) using coordinates <fi n and z n , 



S n = S(^/l - z 2 cos0 n , \J\ - ,z 2 sin <?!>„, z n ) , (4) 
the Hamiltonian becomes 
JS 2 



H = -— £ [v/(l-4)(l-4Jcos( 



1 ) + Xz 



(5) 



where we use a simplified notation for the nearest neigh- 
bors, m = n + a. At low temperature, it is reasonable to 
expand for small z n and (0 n — 4>m)i keeping up to quartic 
order in these, obtaining, 

H = H + H 1 (6) 

where the harmonic and quartic Hamiltonians are 

JS 2 



n 



1 + Az n z n 



(7) 
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(8) 



It is convenient to carry out the calculations in 
wavevector space by defining Fourier transforms 



" q ; 



(9) 



together with a similar definition for Zq. Here q denotes 
a two-dimensional wavevector, the sums are over the first 
Brillouin zone (BZ), and N is the total number of lattice 
sites. Then we rewrite Hq and Hi as 

H = 2JS 2 K 1 " 7q)0q^-q + (1 ~ A 7 q)z q 2f_ q ] , 

q 

(10) 



and 
Hi 



JS 2 

N 



E {50-*- 



^ +q 2 y z qi^q2 Z q3 Z q4 

qi,q2,q3,q4 

(1 + 7q 3 +q 4 - 7q 3 - 7q 4 ) 2 qi Z q 2 0q 3 0q4 
1 

gl 1 + 7qi+q 2 + 7qi+q 3 + 7qi+q 4 



7qi - 7q 2 - 7q 3 - 7q 4 j^qi <P«Ja 0q 3 <Pq4 



} ■ (11) 



For the square lattice, with coordination number z = 4, 
the 7 q function 



7q 



= ~^V qa = ^(cOSfe +COS9 y ) 



(12) 



is an even function of q, i.e., 7 q = 7 -q. The sums in (|l 
are constrained to qi + q2 + q3 + q4 = 0. 

Using the following canonical transformation to cre- 
ation/annihilation operators, 

4>q = a q (a q + o-q), Zq = «/3 q (a q - a_ q ) , (13) 



1 - ^7q 



4S 2 (l-7q) 



1/4 



/?q = 



l-7q 



45 2 (l~A 7 q) 



the quadratic Hamiltonian is diagonalized, 



1/4 



(14) 



(15) 



w q = 4J5^(l- 7q )(l-A 7q ). (16) 
The diagonal parts of Hi can then be written as 



Hi 



^(qi>q2,q3,q4) a qi a q2 a _q 3 a_ q4 , 

qi+q 2 +q 3 +q4=0 



(17) 



where 



A (qi,q 2 ,q3,q4) 



JS 2 
N 



/5 qi /3q 2 /3q3/3q 4 (3 7qi+q 2 7qi+q 3 7qi+q4) 

— 3/3 qi OCq 2 a q3 /3q 4 ( 1 + 7q 2 +q 3 — 7q 2 ~ 7q 3 ) 

— /3q 2 aqiaq 4 ^q 3 (l + 7qi+q 4 - 7qi - 7q 4 ) 
+ /3 q3 a qi CXq 2 /3q 4 ( 1 + 7qi+q 2 — 7qi — 7q 2 ) 
+ Pq 1 0iq 3 CXq 4 f3q 2 (1 + 7 q 3 +q4 — 7q 3 ™ 7q4 ) 

— aq 1 Q!q 2 aq3aq 1 (l + 7qi+q 2 "F 7qi+q 3 F 7qi+q4 

~ 7qi - 7q 2 - 7q 3 ~ 7q 4 ) } ■ (18) 

We have dropped terms in Hi whose expectation values 
in the ground state are zero. Note also that if a par- 
ticular qi is zero, then a qi diverges, however, when this 
occurs the corresponding factors multiplying a qi go to 
zero faster, canceling the divergence. 

III. SPINWAVE DAMPING DUE TO Hi 

The time evolution of the creation/annihilation oper- 
ators as given from the harmonic Hamiltonian is 



a q (t) = a q (0)e l 



aq{t) = a q (0)e~ 



(19) 



As discussed in the Introduction, the effect of the Hi 
perturbation is to shift the spinwave frequencies, in the 
most general case, by some complex self-energy S q . For 
the moment, we ignore the real part of S q , whose effect 
will be included by the memory function approach, and 
concentrate instead on the contributions to the imagi- 
nary part of S q , which determine the spinwave damping. 
To find E q , the contributions of Hi to the single particle 
Green's function will be evaluated, keeping only the di- 
agrams that have a nonzero imaginary part. The lowest 
order diagrams that satisfy this requirement are found to 
be second order in Hi ; there is no contribution from first 
order diagrams. _ _ 

Following the standard procedureJj~u the time-ordered 
zero-temperature Green's function, in second order, is 



iG 2 (k 2 ,ki,i2 - ti) 



H) 2 

2! 



2 roc 



dt'i I dt' 2 X 



(^olTiHiit^H^aMaliti)}]^ . (20) 

Here T is the time ordering operator, and $o is the 
ground state of Hq. Physically, this Green function rep- 
resents the amplitude for a quasiparticle (i.e., spinwave) 
created at time t± in state ki to be destroyed at time t- 2 
in state k 2 , after undergoing two scattering events due 
to Hi. The free-particle Green function, 
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iG (k 2 ,k 1) t 2 -t 1 ) = ($o|T{a k2 (t 2 )4 1 (*i)}l $ o) , (21) 

involves no scattering. Due to momentum conservation, 
these Green functions are proportional to the Kronecker 
delta, <5ki,k 2 j so we can drop the double momentum in- 
dices. Also, because of time translational invariance, they 
depend only on the time difference t = t 2 — 1\. Introduc- 
ing the time Fourier transforms, 



G n (k,t) = ±- / duj G„(k, 



(22) 



the free propagator becomes the well-known expression 

1 



G (k,w) = 



uj — wk + iS ' 



(23) 



where S — + . 

We apply the rules for Feynmann diagrams to evaluate 
the full propagator, G(k, uj), in terms of Go and G 2 . The 
free propagator Go (k, uj) is represented as a straight line 
with an associated momentum k and frequency uj. The 
interaction term in ( |l7| ) is represented as a wavy line as 
in Fig. |l|, with an outgoing straight line (free propagator) 
for each creation operator and an incoming straight line 
for each destruction operator, both with assigned mo- 
mentum and frequency. The two leading diagrams that 
contribute to G 2 (k, uj) are shown in Fig. g. Following the 
standard Feynmann rules for these diagrams, with sum- 
mations over the intermediate frequencies and wavevec- 
tors, and using momentum and frequency conservation 
at each junction, G 2 (k, uj) can be written as 



G 2 (k,w) = G (k,^)S(k, w )Go(k, w ) , 



(24) 



where the self-energy function due to the scattering parts 
of the diagrams is 



E(k,a,) 



duj' f duj" 



The location of the pole gives the perturbed quasiparticle 
energy, u>' k = uj-^ + Re£(k, w k ), and the damping = 
r(k, w k ) is given by 



T k = -Im E(k,w k ) 



(28) 



To apply the result to the system at finite temperature 
T, and obtain the self-energy associated with the temper- 
ature Green's function, the integrals over frequency are 
transformed to sums according to 



^*V) - If; f^) , 



(29) 



2ir 1 

p k B T 

Thus the frequency integrals in (^) become 
du)> f du," 

7; — / -7; — <jo ^0 (jo 

Z7T J Z7T 

_ [n(uJ q ) - re(-^ k+p _ q )][n(a; q + a; k+p _ q ) - n(uj p )} 

W - (^k+p-q +W q - Wp) 

where the Bose factors are 

n(w) = (e*" - l)" 1 . (32) 

The latter form applies to the classical limit. The imag- 
inary part of X is derived from the analytic relation, 



(31) 



1 



x + iS 



X 



inS(x) 



(33) 



G (q, uj") G (p, uj' - uj) G (k + p - q, uj' - uj") 
A(q, k + p q, -k, -p) x 

[^(k,p,-q,-k-p + q)+ J 4(p,k,-q,-k-p + q)] . (25) 



where P takes the principal part and 5(x) is the usual 
delta-function. Because q and p each have two compo- 
nents, the sum in ( p5| ) is equivalent to a four-dimensional 
integration, once we take the continuum limit. Each di- 
mension converts as, for example, — > J dq x , 
where L x L y — Na 2 is the system size, in terms of the 
lattice constant a, leading to 

N 2 a 2 



EE 

q p 



(2tt)4 



dq x \ dq y 



dp x 



d Py . (34) 



S(k, cj) is a sum over the "direct" and the "exchanged" 
diagrams; the factor of 4 is the multiplicity of each of 
these under rearranging the exterior lines. S(k, uj) is seen 
to be G 2 (k, uj) with its incoming and outgoing lines re- 
moved. 

The full interacting propagator G(k, uj) can be approx- 
imated as a sum over repeated applications of SGq, i.e., 



The factors of TV here are canceled by those in the expres- 
sion for A(qi, q 2 , q3, q4), leading to a size-independent 
result. Thus there finally results for the damping func- 
tion, 

T(k, uj) = J d 2 q J d 2 p F(q, p) 6(w - ^k+ P - q + ^ P - w q ) 

(35) 



G(k,w) = Go + Go(£Go) + G (£G ) 2 + .. 



(26) where the kernel F(q, p) is given by 



which, after summing, leads to the well-known Dyson 
equation, 



G(k,w) 



1 



1 



Go -£ uj-uj k + i5- £(k,u>) 



(27) 



F(q,p)=47iY^ A(q,k + p-q,-k,-p) x 

[^(P, k, -q, -k - p + q) + A(k, p, -q, -k - p + q)] x 
[n(w q ) - n(-cj k +p-q)][n(w q + w k+p - q ) - n(w p )] . (36) 
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We use this result to calculate the damping at the spin- 
wave frequency: Tk = r(k, u> = Wk), in the thermody- 
namic limit L — > oo. Note that in the classical limit, 
the Bose factors are simply proportional to temperature. 
Therefore, we only need to calculate I\ at one tempera- 
ture (T = 1.0 J) and then scale with T 2 to get its value 
at other temperatures. 



Shifting coordinates to u = x — x c , and writing the gra- 
dient at the cell center as Vg(x c ) = |Vg(x c )|l, where 
the unit vector 1 = (h,h, h, h) points along the gradient 
evaluated at the cell center, we have 

DceiiM = JpArt I d'u 8(w-u-\), (41) 
|Vs(Xc)| J cell 



A. Evaluation of contour integrals in four dimensions 

The summation/integration in ([55]) involving a delta 
function presents an interesting mathematical problem in 
its evaluation. General integrals of this form (i.e., with 
one delta-function) in d dimensions can be considered, 



D{u) = / d d x f(x)6(u - g(x)), 



(37) 



where x is a d-dimensional vector. In two dimensions, for 
example, the delta-function will pick off a particular con- 
tour that contributes to the integral; in three dimensions 
that contour will actually be a surface, and in four dimen- 
sions it will be a three-dimensional volume. Depending 
on the form of the function g(x), the contributing region 
need not be simply connected, but may consist of sep- 
arated distinct parts. For the higher dimensional cases, 
it is difficult to visualize the region or regions that con- 
tribute to the integral. Thus it is necessary to have a 
way to evaluate the integral without appealing directly 
to geometrical properties. _ _ . 

In a series of papers, Gilat and co-workersBilJ have 
considered how to evaluate these integrals in three di- 
mensions, in work related the calculation of spectral den- 
sity functions. We have used the "linear analytic" (LA) 
methodE3 of Gilat and Raubenheimer, extended to four 
dimensions, which gives good precision for few sample 
points, at little cost in complexity. Suppose we calculate 
D(u>) over a small range of u> and then take the limit 
w -> Wk- To do so, a 4D cubic grid of small cells of size 
(2b) 4 is set in the x-space. Then we consider the con- 
tribution of each cell to D(u>), thereafter summing over 
the cells. Suppose the cells are small enough so that the 
kernel /(x) has only a small variation within a cell, and, 
therefore, we can use its value at the cell center, x c . Then 
we can write 



D(u) = J2D cen (u) 



cells 



Dceu(w) = /(x c ) / drx 5(u> - g(x)) 

J ceil 



(38) 



(39) 



The idea of the LA method is to expand g(x) around the 
cell center, 



g{x) w g(x c ) + (x - x c ) • V#(x c ) 



(40) 



g - fl( x c) 

|Vff(x c )| 



(42) 



The variable w is a rescaling of w. If we use the integral 
representation of the delta function, in the form, 

5(w-u-l) = — dk e **(«-«iii-»afa-«8is-«.i4) ( 



(43) 



then the contribution from one cell is found to be 

4 ,b 



-Dcell(w) 



/(x c 



1 



dk e lkw 



|V S (x c )|2^ J_ 

/(x c ) 1 

|V<?(x c )| I2\hhhh\ £j 



n 

3=1 



-b 



da, e ikU]l] 



16 



j2(-iy m \w+ Wl \ 3 , (44) 



where the Wi are the values of parameter w at the 16 
corners of the cell, i.e., 



b(±h ±l 2 ±l 3 ± U) , 



(45) 



with all 16 possible choices of the signs, and m$ is the 
number of minus signs used in each Wi . One can note that 
similar simple formulas result for this type of integral in 
other dimensions as well .113 

In the actual application of Eq. (Q), inside the cell, the 
parameter w will range between minimum and maximum 
values: 



= H\h\ + \l 2 \ + \h\ + \l4\) = -i 



(46) 



These values hold at two opposite corners of the cell. By 
Eq. (^2|), there will be a corresponding range of ui over 
which this cell contributes to D(u). Then for the general 
problem where D is desired over a large range of ui, what 
is done is to set a uniform grid of w-values, finding the 
contribution of each cell to D(u>) for each value on the 
grid, say, lo = u> n . Typically it works well to have the 
grid fine enough so that each cell contributes to D at 
many different oj n . 



B. Numerical calculation of I\ 



For the damping problem, Eq. (p5[), we have the re- 
placements, 

/(x) = F(q, p) , g(x) = cj k+p _ q - uj p + w q , (47) 
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with the 4D variable x = (q, p). We want to find r(k, us) 
right at lu = lu^ for some particular wavevector, k. To 
accomplish this, we set a range w m i n < oj < w max , and de- 
fine the grid of values uj n = uj min + ra(w max - Lu min )/N gridl 
using A^grid = 500. For the integrations over q and 
p, the cell sizes actually can be quite large: we used 
the cell width 2b = 2ir/L, where L = 41 is a "sys- 
tem size," which gave results consistent with those us- 
ing much smaller cells (larger L). That is, each axis, 
q x , q y , p x , p yi was partitioned into 41 cells. First, to 
get good precision, a narrow range symmetrical around 
the spinwave frequency was used: cJ m i n = uj-k — 0.1JS, 
w mM = c^k + Q.US. Second, for comparison, the behav- 
ior of the damping function was calculated over the whole 
range of possible spinwave frequencies, i.e., w m i n = 0, 
w max ~ cj q= ( Wi7r ) = 4-^/2(1 + X)JS. In both of these 
approaches, interpolation on the u-grid gives r(k, u>k). 
The scheme was found to be very robust, with both ap- 
proaches agreeing to better than one percent precision. 

We should also note one other important calculational 
detail. Because Eq. (^) depends on the reciprocal gra- 
dient of g, it is necessary to avoid any points in the q, p- 
space where this gradient vanishes. We accomplished 
this by taking wavevectors k = (xk,Uk)^-, (at L = 41), 
and then shifting the q, p-cell center positions slightly 
from the usual values: q = (x q + \,y q + 5)^7 + (^j A), 
P = (x p + \,y p + \)^f + (A, A). Here x k ,y k ,x q ,y q , 
etc., are integers ranging from to L — 1, and we used 
A = 0.022378 to shift the q, p-cells away from high sym- 
metry points. The choice of an odd number for L also 
helps to avoid certain high symmetry points. Some slight 
improvement in precision can be achieved with a larger 
value of L, however, at the expense of a tremendous in- 
crease in cpu time, which varies as 4 L . At L = 41, the 
typical relative errors in r(k, lu) were found to be much 
less than one percent, when compared to calculations us- 
ing L = 81. 

Some results for the damping function, which is pro- 
portional to T 2 , are shown in Fig. [|for A = 0. In general, 
it is found that T(k,a;) either has a cusp or a peak at 
the spinwave frequency (arrows in Fig. ^). The damping 
strength increases with |k|, as seen in Fig. [|. Along the 
(11) direction, however, Tk tends to saturate halfway out 
in the zone. Keeping in mind that Fk oc T 2 , the damping 
strength we find here is quite weak, especially when we 
consider its effects at T <C J where the spinwave approx- 
imation is valid. 

At small nonzero values of the anisotropy parameter 
A, there will be only a very slight change in Tk, because 
A does not appear explicitly in Hi, but only enters the 
calculation via the dispersion relation and via the a and (3 
functions of (|l3|). To demonstrate this effect, Fig. [3] shows 
Tk at various values of A for k along the (10) direction. 
In general, the damping decreases with increasing A, i.e. 
as the system becomes more isotropic. 

The calculation of the memory function requires the 
knowledge of Tk over the entire BZ. We can reduce the 



region of calculation to k x > k y > 0, because, by sym- 
metry, T k = r_ k , and T (k ^ ky) = T( kv , kx y In this region 
Tk was calculated on a uniform square grid with spacing 
Ak x = Ak y = 27r/50. Linear interpolation was used to 
get values between these grid points. 



IV. THE MEMORY FUNCTION FORMALISM 
WITH DAMPING 

For the calculation of the dynamic correlation function 
S(q,uj), we use the memory function formalism due to 
Reiter and Sj Slander ,EJ who applied the procedure to the 
isotropic Heisenberg model. The technique is based ©a 
a low-temperature expansion of the memory function.Ej 
The formalism was applied to the ID isotropic .ajjitiferro- 
magnet in an external field by De Raedt et alll3 and to. 
a ID easy-plane antiferromagnet by Gouvea and Pires.EEl 
More recently, it was also apalied to 2D easy-plane fer- 
romagnets by Menezes et alJiH although some numerical 
difficulties inhibited numerical evaluation of the analytic 
formulas. Some of thf se problems were surpassed in the 
work of Pereira et al.EXLj In all of these works, no mode 
damping was included. The presence of damping actu- 
ally relieves the numerical difficulty of the problem, as 
we see below. 

The dynamic correlation function is defined in terms 
of the Fourier transformed spin-fields as 

1 r°° 

S aa (q,iu) = —J dt e— t (5 a (q,t)5«(-q,0)) , (48) 

where a is any component of S, and by symmetry, 
S xx (q,uj) = S vv (q,uj). In this and following formulas 
for correlation functions, we note the distinction between 
continuum (functions of q) and discrete (subscripts q) 
Fourier transforms, rWjhich are related by factors of 2n/L 
for each dimensionEJ This is necessary so that we can 
make a careful comparison between the theory and the 
MC-SD simulations, including the absolute magnitudes. 
Specifically, the continuum correlation function in Jis}) , 
which gives a thermodynamic weight over a wavevector 
range d 2 q = (2ir/L) 2 , is related to the discrete corre- 
lation at a specific wavevector «, using discrete Fourier 
transforms as defined in ([)]), bycil 

(£«(q,t)S a (-q s 0)) = -l-<^(t)5« q (0)) . (49) 

In the analytic calculations, many intermediate quanti- 
ties are discrete functions, however, for the comparison 
with the MC-SD data, both results will be converted and 
displayed as continuum quantities. 

Due to the symmetry in the xy-plane, the calculation 
is carried out by making averages of the rotationally 
symmetric quantity, = S x + S%, with S xx (c[,uj) = 
hS (q, u>). Then in the memory function calculation, 
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the continuum dynamical structure function can be writ- 
ten in terms of a second-order memory function, (lu), 
to be described below, according to 



S X - L (q,w) = 

(2 7 r)- 2 7T- 1 (L^,L^ q ) Im S^(fa,) 

p - + w Re E^H] 2 + [ w Im E£(u,)] : 



(50) 



where L = —i d/dt is the Liouville operator, and the 
second moment (w 2 ) q is determined from a ratio, 



(51) 



The location of the spinwave peak that results from Eq. 
( |5C)| ) is very sensitive to the value of (t^ 2 )q , therefore we 
are especially concerned to get an accurate evaluation of 
it. The denominator of (|5l]) is a static correlation func- 
tion, whose value will be taken from the Monte Carlo 
simulations, although for low temperatures it can be eval- 
uated approximately as 



(^» I7 (l-7 q )- 1 



(4) 



T 

4J ' 



(52a) 



(52b) 



The numerator of ( pl| ) has been shownllZHIa to have a lead- 
ing low-temperature dependence, also related to static 
nearest neighbor correlations, 

(LSq , LS_^) = 

4JT(1 - A 7q )(S^) - 8JT( 7q - X)(S^SS) ■ (53) 

The correlation {S^-Sq) oc AT, therefore its contribution 
can be ignored for small A. There being two bonds per 
lattice site, the in-plane near-neighbor correlation is re- 
lated to the average energy per bond by 

( 2Jn) = S * ~ ^ ■ ^ 0> " ^ " {S ^ S ° ] ■ (M) 

This correlation can also be obtained from the MC sim- 
ulations, although a good low-T approximation is 
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— \ e -l-273*7 



(55) 



A comparison of the analytic and MC static parameters 
will be given below. 

To the leading order in temperature, the time- 
dependent memory function has been shown to be given 
from 



E£(t) = -{LS^LS^)- L M^(t) , 



M°(t) = (QL 2 S^e- lLt QL 2 S^ 



(56a) 



(56b) 



where Q is projection onto the non-secular variables, see 
Ref. O, whose action is equivalent to 



QL 2 S« = (L 2 - (loX)S« 



(57) 



Using expressions (|j) expanded to second order in tf> n , 
z n , and (|l3|), it is straightforward to write the spins 
S n in terms of the creation/annihilation operators. In 
the leading order in temperature, the resulting expres- 
sions involve expectation values of products of four cre- 
ation/annihilation operators. The damping is included 
in the time evolution according to, 



't(t) = o t(0)e<* < '- r «>* , a q (i) = a q (0)e(" 



-r,)* 



(58) 



Then, there results, after some manipulations, expres- 
sions to give the memory function, 



M: 



"(*) = W / +(q,P)e- r+(q ^ p)t cosO + (q,p)t 
p 

+IU_(q,p)e- r - (q < p)t cosr2_(q,p)£j . (59) 



The two terms correspond to spinwave sum (f2 + ) and 
difference (f2_) processes, with weights 



W±(q,p) = ii(^ +p )n(wa_ p )(sT! ) 2 , 

where 

S = S/3. +p /3,_ p {^ 2 -(4JS) 2 x 

(l-7H +P )(7a_p-A7q +p ) 
+ ( 1 -7H-p)(7s +p 



(60) 



A7C 



} , (61a) 



t = S a i+p a^ p ^u; 2 + (4JS) 2 x 



(l-A7|_ p )(7|_ p -A7^ +p ) 
+(1- A7| +p )( 7 | +p -A 7 |_ p ) }. (61b) 
The frequencies and damping rates of these processes are 
^±(q,p) =^ +p iw,_ p , (62a) 



r±(q, P ) = ra_ 



(62b) 



Finally, the frequency-dependent memory function is 
given by the one-sided Fourier-Laplace transform, 



/>OC 

Jo 



(63) 



From this expression, Eq. (^), and equations following, 
we obtain a sum over all the sum/difference processes, 
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P 



+W- 



ft. 



iT 4 



uj + ft_ — ir_ 



( w + ft_) 2 + r 2 _ 



U) 


- ft+ 


-zr + i 


(CJ 


- ft + 




U) 


- ft_ 


-ir_ ] 


(UJ 


- ft_ 


) 2 + r 2 J 



}. (64) 



where the (q, p) dependencies of VF± , ft± , and r± have 
been suppressed. This is the expression used to evaluate 
and thus the dynamic structure function, via Eq. 

To evaluate (Q) on a discrete L x L lattice, the 
wavevector q must be one of the allowed values, 



(x q ,y q ){2n/L) 



(65) 



where x q and y q are integers. It is also necessary to make 
the summations over p such that the vectors § ± P are 
both given by integers times 2it/L. Thus, it is convenient 
to change from the summation over p to a summation 
over a new variable k defined by 



k = 



(x k ,y k )(2ir/L), 



(66) 



we also have ^ 



where Xk and y k are integers. Then, 
p = q — k equal to an allowed wavevector. This insures 
that the sums will be over spinwaves that are allowed 
on the given finite lattice. This approach is essential 
for calculation of the finite size dynamical effects. On 
the other hand, for the infinite sized system, no variable 
change is needed, and there is the usual replacement, 

Note that in the limit of zero damping, Eq. (Q) is 
inapplicable, except for the continuum limit (i.e., L — ► 
oo). The resulting continuum integrals for the infinite 
system are singular, and the analytic relation (|3^) must 
be applied, leading to 



I P I d 2 p]\ 



1 



P / d 2 p 



uj + ft + uj — ft_ 
1 1 



ft_ 



- ft_ 



in / d 2 p W+ 6(uj + ft + ) + S(uj - ft+) 



d 2 p 



8{uj + ft_) + 6(uj - n_) 



(67) 



Indeed, in the very low temperature range (T < 0.2J) 
where the damping is extremely weak, it may be more ef- 
ficient to use this expression rather than Eq. ( |64| ) , which 
has poorer convergence with L. Eq. (j6?j) was applied by 
Pereira et alJi§E3; the imaginary part can be evaluated by 
techniques described here (Sec. ill A), and the princiftal- 
valued real part can be calculated by an extensionES of 



techniques used in ID £3 In Ref. [19] the last expression 
was found to give finite amplitude spinwave peaks, with- 
out damping present. 

Before presenting the results obtained via the memory 
function approach, the key aspects of the Monte Carlo 
and spin-dynamics simulations are summarized. 



V. MONTE CARLO AND SPIN-DYNAMICS 

To test the validity of the memory function technique, 
we compare it with numerical simulations on L x L square 
lattices (L — 128, A = 0) using Monte-Carlo and spin- 
dynamics (SD) simulations, which include effects due 
to all thermodynamically allowed excitations. As men- 
tioned above, we also use some data for static correla- 
tions from the Monte Carlo calculations as input to the 
memory function dynamical calculations. The techniques 
used here have been described in Ref. and are based 
on the simulation methods developed in Ref. 2^. 

In the MC calculation, we applied a combination of 
Metropolis single-spin moves and over-relaxation moves 
that modify all three spin coraoonents, and in addition, 
Wolff single-cluster operationsE3 that modify only the xy 
spin components. The over-relaxation and cluster moves 
are important at low temperatures, where the xy spin 
components tend to freeze and single spin moves become 
inefficient. 

In the combined MC-SD simulations, a set of MC gen- 
erated states from thermodynamic equilibrium at tem- 
perature T are used as initial conditions for numeri- 
cal integrations in time via a fourth order Runge-Kutta 
scheme. Specifically, the first 4000 MC steps were used 
for equilibrating the system. Then, every 400 MC steps 
a state was taken to initiate the SD integration. This 
was repeated 500 times; the displayed 5(q, uj) data are 
averages over the 500 time integrations. The total MC 
sequence was then 2 x 10 5 steps; every 25 steps a state 
was used to calculate averages of the static correlation 
functions. 

In the time integration, the basic Runge-Kutta time 
step was At = 0.03( JS)" 1 , and every 11 time steps, 
data samples (for Sq(i), etc.) were saved for q's in 
the (10), (01), and (11) directions only (due to ma- 
chine memory restrictions). A total of 2 12 data sam- 
ples were saved, followed by applying a fast-Fourier- 
transform to obtain S xx (q,u>). The total time of inte- 
gration of i max = 1351.68(JS') _1 is large enough so that 
no smoothing window is necessary for the data genera- 
tion. Then the frequency resolution of these results is 
good, Auj = 2n/t max = 0.0046484 JS, so that any fine 
details in the low temperature results should be visi- 
ble. For the lowest temperatures (T < 0.3 J), where 
we wanted to see the complex fine details due to the fi- 
nite size effects, some simulations were made saving data 
samples every 28 times steps, leading to a total integra- 
tion time L a 



3440.64(JS f ) , and frequency resolu- 
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tion Aiv = 2vr/i max = 0.0018262JS'. 

VI. NUMERICAL RESULTS 

First, a comparison of the low-T formulas and static 
MC results needed as input to the memory function for- 
malism is presented. Then dynamical results from the 
memory function and from MD-SD will be presented for 
large systems, and finally, for smaller systems, where the 
finite size features will be analyzed. 

A. Statics 

The continuum static structure function, 

S* x (q) = (S^Sf q ) , (68) 

was measured directly in the MC simulations. A com- 
parison with the results based on the low-temperature 
expression (62) is given in Fig. | There we also show 
a fit of the MC data to a form similar to the analytic 
expression, but with a fitting coefficient (A) and power 
(B), 

(Sq S- q )fittcd = A(l - 7q r B • (69) 

From the fits, the parameter A is somewhat larger than 
the low-temperature limit, A = T/(4J). This can be ex- 
pected due to the renormalization of the effective J to 
lower values with temperature. The power B falls with 
increasing T, which must also be due to anharmonic ef- 
fects. There is little difference between the low-T expres- 
sion and the MC data at the lowest temperature. The 
deviations between the MC and low-T theory increase 
substantially by T = 0.5 J, especially at higher wavevec- 
tors. Note that the fits are only shown for demonstration; 
we used the MC data itself as inputs to the memory func- 
tion calculations. 

The other static input to the memory function calcu- 
lation, for determining the second moment, is the quan- 
tity (LSqLS^), which is related to the nearest neighbor 
correlation (S^-Sq) by Eq. (f53]). Estimates of that latter 
from the MC evaluation of the average energy per spin 
[Eq. ( |54[) ] are compared to the theoretical expression ( p5| ) 
in Fig7[7|. The analytical expression is seen to give results 
consistent with the MC data up to T w 0.5 J. 

When combined to give the second moment, the MC 
data produces results as shown in Fig. |8[ It is essentially 
this quantity that has the strongest influence on the lo- 
cations of the spinwave peaks in 5(q, u). These are the 
data used as input to the memory function calculation. 



B. Dynamics: Memory Function compared with 
MC-SD 

We evaluate Eq(w) by using the sum over discrete 
wavevectors p as displayed in Eq. (J64|), and then apply 
Eq. © to determine S xx {q,uj) = jS ±± (q,Lj). If L is 
taken large enough, then the sum should be a good ap- 
proximation to the infinite size limit. In general, we can 
see whether larger L is necessary by the presence of finite 
size features in the data for Eq and S xx (q, ui). Typically 
we used L = 1024 to give a reasonable approximation to 
the continuum limit for T > 0.3J and larger q-values; at 
lower temperatures or wavevectors even larger L is nec- 
essary to smooth out all the finite size features, due to 
the very weak damping. 

In Fig. [9| some typical results for the real and imagi- 
nary parts of Eq (w) are shown for A = 0, T — 0.3 J, and 
q = (tt/4, tt/4), where L = 2048 was used in the sum. 
The resulting dynamical structure function is compared 
to the corresponding MC-SD data for L = 128 in Fig. [H]. 
At the (zero-T) spinwave frequency u> q = 2. 16478 JS, Im 
Eq (a;) achieves its minimum value, which is not zero due 
to the presence of the damping, while Re Eq (u>) has a 
small nonzero value there as well. The spinwave peak in 
the memory data occurs at u) = w pC ak = 2.005 JS, which 
is just slightly greater than the associated value of the 

second moment, lu^ = J (w 2 )^ = 1.9648JS'. The spin- 
wave peak position will occur where the denominator of 
Eq. ( |50| ) achieves its minimum value, which gives, to a 
good approximation, 

^peak = < ^/l-ReEq-«)/< . (70) 

For this case, the memory function data give 
ReEq(o>q) w —0.081JS, which is consistent with this 
relation. The peak position is very sensitive to the value 
of the second moment, which is why we have tried to esti- 
mate it as accurately as possible from the static MC data. 
Still, the MC-SD data show the spinwave peak closer to 
uj = 1.98JS, slightly lower than the memory function 
prediction. Using Eq. (|50|), the spinwave peak height is 
proportional to [ImE(w pca i()] _1 , while the peak width is 
approximately equal to ImE(c i ;p 0a ] 5 ). When compared to 
the MC-SD data, the memory function calculation usu- 
ally overestimates the peak height and underestimates 
the peak width, although both are in reasonable agree- 
ment with the MC-SD data. In fact, these errors tend to 
be larger at lower temperature. 

In Figs. p~T| — |l3|, we display some further comparisons of 
the memory function calculations and MC-SD results at 
other values of q for temperatures ranging from T = 0.1J 
to T = 0.5J. For the lower temperatures there is very 
good agreement between the memory theory and MC- 
SD peak positions, while the peak heights and widths do 
not agree as well. On the other hand, at higher tem- 
peratures the opposite is true: the memory data show 
greater errors in the peak positions, especially at large 







g, but the peak heights and widths are predicted very 
well. Probably these discrepancies are caused by a fairly 
large underestimate of the damping strength at low tem- 
perature, as well as errors in the estimate of the second 
moment at higher temperatures. However, in general, 
the tails of the spinwave peaks far from w q are described 
very well by the memory function calculations. 

C. Dynamics: Small Systems and Finite Size Effects 

When we consider the results using relatively smaller 
values of wavevector, S xx (q,uj) shows not only a spin- 
wave peak, but additional features (small peaks) which 
depend on the system size. These features are obvi- 
ous in the MC-SD S(q, u/) data shown above. Similar 
size-dependent dynamical features have bee n notic ed in 
earlier MC-SD simulations of XY-symmetryESEZIGa and 
easy-axis symmetryc£l magnets, but without a clear the- 
oretical explanation. We now consider a closer look at 
the details of these features, making a comparison with 
the memory function calculations at smaller values of L. 

Some examples showing the resulting memory func- 
tion data for S xx (q,u>) for L = 128 systems are given 
in Fig. 0a [q = (16,0)2tt/L, T = 0.3.7], Fig. 0a 
[q = (9, 0)2tt/L, T = 0.3 J], and Fig. 0a [q = (5, 0)2tt/L, 
T = 0.1 J]. In fact, there appear various sequences of 
small peaks. The number of separate sequences appears 
to depend on the wavevector in integer units of 2tt/L. 
Note that the logarithmic scale exaggerates the size of 
these subpeaks, they are typically a few orders of mag- 
nitude smaller than the spinwave peak. For sufficiently 
low temperatures, these effects also ap pea r clearly in the 
MC-SD results: see Figs. 0b, 0b and [0b, although the 
corresponding features there are smeared out compared 
to the memory function calculation. The effects are also 
present at larger wavevectors, however, in that case there 
are many more of these finite size subpeaks, which tend 
to merge together into a smoother background. 

These extra peaks can each be associated with partic- 
ular spinwave sum or difference process. To see that this 
is the case, consider first the difference processes, whose 
frequencies Q_(q, p) are always less than the spinwave 
frequency cj q . A plot of the weights W- [Eq. (0)] that 
determine the contributions to the memory function ver- 
sus the corresponding frequencies f2_ = uj^ — w q _k, can 
show how strongly each process [i.e., each k as defined in 

f)] contributes to Im Eq and hence to S(q, uj). In Fig. 
we made such a plot for an even integer q, namely, 
q = (16, 0)(2tt/Z), for L = 128 and T = 0.3J. Each dot 
there corresponds to one allowed wavevector k; some of 
the values of k in units of 2ir/L are indicated there. The 
sequence of k's, (8,0), (8,1), (8,2), etc., leads to multi- 
ple contributions at zero frequency, that is, a weak cen- 
tral peak, which is a general result for any even integer 
q. The sequence, (9,0), (9,1), (9,2), etc., produces a se- 
quence of small peaks in 5(q, uj) out to w »s 0.2JS, and, 



the next sequence, (10,0), (10,1), (10,2), etc., gives a se- 
quence of peaks out to u m 0AJS, and so on. Each 
of these k-sequences can be identified with correspond- 
ing peak-sequences in the memory function calculation of 
5(q, cj), as seen in Fig. 0a. Indeed, for the wavevector 
q = (16, 0)(2n/L), there are 8 + 1 of these sequences be- 
low the spinwave peak, including the ones that contribute 
only at w = 0. The higher sequences are more spread 
out, and merge with the other ones, making a compli- 
cated fine structure in S*(q, uj). Similar features are seen 
in the MC-SD data, Fig. 0b, although it is clear that 
the individual peaks of each sequence are smeared out, 
resulting in one wider feature from the entire sequence. 

Another example of these finite size features, for an 
odd integer wavevector, q = (9, 0)(2w/L) is shown in 
Fig. [U§ again for L = 128 and T = 0.3 J. For this smaller 
integer wavevector there are now only 5 sequences of fine- 
structure difference process peaks. It is significant that 
there is no such peak at zero frequency. This is clearly 
caused by the absence of any allowed spinwave difference 
process exactly at il_ = 0, although there are many such 
processes with very weak weights close to f2_ = 0. In 
general, for any odd integer wavevector, <S(q, uj) will have 
a relative minimum at uj = 0, because of the absence of 
allowed difference processes there. 

In Fig. we also show some of the weights of the sum- 
processes, i7+ = cjk +^q-k, indicated by plus signs. For 
either sum or difference processes, numbers in parenthe- 
sis refer to the values of k in units of 2ir/L. The sum 
frequencies must fall above the spinwave frequency u q . 
Comparing Figures and 0, we see that the finite size 
features above w q in 5(q, uj) can each be identified with 
a particular spinwave sum process. The identification, 
however, is somewhat more complex than that for the 
difference processes, because the different sequences of 
processes overlap in frequency. The match between the 
sum peak frequencies in the memory data and the MC- 
SD data is not as good as for the differences processes 
below ujq. We can note that at these higher frequencies, 
any relative errors appear larger on an absolute scale, 
and, the spinwave sum/difference frequencies do not in- 
clude any exchange softening effects. 

Finally we show another example, for q — (5, 0)(2ir/L) 
at T = 0.1J, L = 128 in Fig. 0, where both the differ- 
ence and sum processes are shown. In this case there are 
only 3 dominant sequences of difference processes leading 
to finite size features. 



D. Higher Temperature 

Although the memory function formalism applied 
above, strictly speaking, applies to low temperatures, it is 
interesting to consider what kind of predictions it makes 
if applied at higher temperatures. It is clear from the 
above comparisons between MC-SD data and the mem- 
ory function data that at low temperatures, where the 
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spinwave approximation should be valid, the damping of 
the modes has been underestimated. This is seen espe- 
cially in the T = 0.1J memory data where all the fi- 
nite size features, as well as the spinwave peak, are too 
narrow. Thus it is interesting to see whether the spin- 
wave width and shape are better described by this theory 
at higher temperature, because the theoretical damping 
rates grow as T 2 . The other reason to look at higher 
temperature is to consider how the background inten- 
sity in 5(q, to) for frequencies u> < w q , depends on the 
temperature. For one thing, the increase in tempera- 
ture should spread all of the finite-size effect peaks and 
smooth the spinwave background, but can be expected to 
raise its intensity as well. For comparison, as T crosses 
the vortex unbiading-, temperature Tbkt, previous MC- 
SD simulationsaO-ES show that this background grows 
with T, until it dominates over the spinwave peak. The 
effect has been described by a theory employing a gas of 
freely moving vortices,E3 although simulations have given 
some indication that vortex lifetimesci are rather too 
short to take this theory literally, and, vortices tend to be 
pinned by the lattice and also have complex motions not 
resembling free translation. Thus it is interesting to com- 
pare the memory approach to the MC-SD data at higher 
temperatures, where the latter will include effects due to 
all possible types of excitations, beyond linear spinwaves. 

For A = 0, the vortex unbinding temperature has been 
estimated to be Tbkt ~ 0.7J. In Fig. gfj, we show some 
comparisons of the continuum (L > 1024) memory and 
MC-SD [L — 128) data for S(q, u>), at temperatures 
0.5J > T < 0.9J for two particular values of q. As 
already seen, the spinwave peak locations in the MC-SD 
data here also fall lower than the memory function theory 
prediction. Also, there is a significant growth of the low- 
frequency background (u> < w q ) in the MC-SD data and 
in the memory function calculations. For T = 0.9J the 
low-frequency background already dominates over the 
spinwave peak, if the latter can be seen at all. The effect 
is stronger in the MC-SD data: there is a clear excess 
of intensity around zero-frequency in 5(q, uj) in the MD- 
SD data compared to the memory function data. This 
effect is absent at T < 0.5J, however, where the MC-SD 
and memory data coincide for low frequencies; in this low 
temperature region the spinwave approximations used in 
the memory approach are valid. For temperatures near 
and above Tbkt, this extra central peak intensity, not de- 
scribed within the spinwave picture of the memory func- 
tion approach, may possibly be attributed to nonlinear 
excitations, including specifically, vortices. This is con- 
sistent with the phenomenological vortex gas description 
of Mertens et al:E3 the excess central peak component is 
stronger for lower wavevectors, because the vortex pairs 
being generated above Tbkt disorder the spins on long 
length scales. 



VII. DISCUSSION AND CONCLUSIONS 

The presence of spinwave damping included into the 
memory function approach has been suggested here as a 
way in which finite amplitude spinwave peaks result in 
the 2D XY model, canceling the divergences in previous 
theories. A rough calculation of the damping strength 
has been made using the Feynmann diagram technique 
applied to the quartic corrections of the harmonic Hamil- 
tonian. The diagrams used involve only two scattering 
events; clearly this gives an underestimate of the damp- 
ing, and it is likely that including the third order dia- 
grams could make a substantial modification to the re- 
sults. Generally the damping Tk obtained increases away 
from the center of the Brilluoin zone, being close to zero 
near the center and being greatest along the (11) direc- 
tion. In this classical calculation, the damping varies as 
the square of the temperature; this results in extremely 
weak damping strength (Fk <C t^k) for temperatures be- 
low 0.3J. As the temperature approaches and passes 
Tbkt, the damping becomes comparable to the spinwave 
frequency. 

The comparison of the memory function data and MC- 
SD data for S"(q, ui) is good in several aspects: the theory 
gives finite spinwave peaks, they are positioned correctly 
for low temperatures, and the multiple finite-size-effect 
peaks are reproduced as in the MC-SD data. The most 
significant problem in the memory function approach is 
the extremely narrow peak widths-clearly caused by an 
underestimate of the true damping strength, especially at 
lower temperatures and wavevectors. For very low tem- 
peratures in the present theory, the longest wavelength 
spinwaves will not scatter even once before crossing the 
entire system. Of course, this is a quantum picture, ap- 
plied to the classical system, so perhaps a better ap- 
proach in this low-T limit would be to consider some 
kind of adiabatic scattering in the nonuniform potential 
due to the spinwave population. This would essentially 
be the equivalent of the quantum calculation at very high 
order diagrams rather than at the lowest order. 

When the memory calculation is carried out for large 
systems, that is, in the continuum limit, the finite-size 
peaks merge into a smooth background. Indeed, in this 
limit, the spinwave peak itself has contributions on the 
high frequency side from spinwave sum processes, and on 
the low frequency side from spinwave difference processes. 
The combination of all of these contributions, essentially 
an infinite number of sum and difference processes, leads 
to an overall smooth and nondivergent spinwave peak. 
Even in the limit of zero damping, the spinwave peak are 
nondivergent, according to-eontinuum calculations [using 
Eqs. ^37j by Pereira et al.lia However, when damping is 
included, the heights of the spinwave peaks are smaller, 
and we expect the results to be more realistic. 

In the memory function theory, the spinwave peak 
height was seen to be proportional to [ImS(wp Ca |()] _1 . 
The function S(q, to), to a leading approximation, is pro- 
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portional to the temperature T. From Eq. (pOj) , in regions 
away from the spinwave peak S(q, uj) should be propor- 
tional to X(q, uj). Thus we see that the memory func- 
tion theory demonstrates how the background grows with 
temperature while the spinwave peak relatively tends to 
fall. For high enough temperature, the distribution of 
spinwave difference processes near zero frequency domi- 
nates the spectrum of S(q,uj), leading to a central peak 
stronger than the spinwave peak. Clearly this means ap- 
plying the theory well above the limiting temperature 
where we expect spinwave approximations to be valid. 
However, this gives results very similar to those in the 
MC-SD simulations, although the latter can be expected 
to include nonlinear effects such as vortices. In fact, we 
have seen a strong indication for nonlinear effects (be- 
yond spinwave difference processes) above Tbkt as extra 
intensity in S^q, uj) near zero frequency in the MC-SD 
data. 

In conclusion, application of the memory function for- 
malism, including a leading approximation to the mode 
damping, has been shown to accomplish three things for 
the XY model: i) produce nondivergent spinwave peaks; 
ii) explain the finite size features in S*(q, uj) both above 
and below u q ; iii) explain the background in 5(q, uj) 
below w q for the infinite system. Note that strictly 
speaking, damping is necessary only for item ii) , because 
5(q, uj) cannot be calculated in the memory function for- 
malism for a finite system without damping. The con- 
tinuum memory function approach gives finite size spin- 
wave peaks and a low-frequency background even in the 
absence of damping. The theory also suggests that a siz- 
able contribution to the growth of the central peak for 
T above Tbkt may be due to spinwave difference pro- 
cesses. The excess central peak intensity in the MC-SD 
data, when compared to the memory function data, must 
be due to nonlinear excitations, and gives evidence for 
dynamical effects due to vortices. 
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FIG. 1. Notation for the representation of the interaction 
terms in Hi, Eq. (jfj]), in the Feynmann diagrams. 
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FIG. 2. The two leading diagrams that contribute to 
G 2 (k,w), as in Eqs. (|3) and (El). 
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FIG. 3. Typical numerical results for the damping func- 
tion r(k,cj) at A = 0, at various wavevectors as indicated. 
The arrows mark the positions of the spinwave frequency Uk, 
for each associated wavevector. 
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FIG. 4. The damping rate T k for A = in the (10) and 
(11) directions. 
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FIG. 6. Static correlation function for A = 0, calculated 
from MC simulations for L = 128 (symbols), with fits to 
formula (pSj) (solid curves), compared to the low-temperature 
expression ( |52| ) (dashed curves). 
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FIG. 7. Static nearest neighbor correlation function for 
A = 0, calculated from MC simulations for L — 128 using Eq. 

(54) (symbols), compared to the low-temperature expression 

( 55) (solid curve). 
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FIG. 8. The second moment vs. q for A = 0, at the 
temperatures indicated, using the data from L = 128 MC 
simulations, compared with the zero-temperature dispersion 
relation (solid curve). 
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FIG. 10. The dynamical structure function obtained from 
the memory function of Fig. | (solid curve), calculated us- 
ing L = 2048 for T = 0.3J, q = (7r/4,7r/4). compared with 
MC-SD data (symbols) for the L = 128 system. The fre- 
quency resolution of the SD data here and in Figs, [ll], [l^, |l3| 
is Alu = 0.0046484JS. 
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FIG. 9. Real and imaginary parts of the memory func- 
tion Eq(w), calculated using L = 2048 for T = 0.3 J, 
q = (7r/4, 7r/4), and the resulting dynamical structure func- 
tion S M (q,w). 
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FIG. 11. The T = 0.1J dynamical structure function 
obtained from the memory function calculations (solid, us- 
ing L = 2048-4096) and MC-SD simulations (dotted) for the 
L — 128 system, at wavevectors indicated in parenthesis. The 
memory data and the MC-SD data are both smeared by fi- 
nite-size effects, even for these large L. 
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FIG. 12. The T = 0.3J dynamical structure function 
obtained from the memory function calculations (solid curves, 
using L = 1024-2048) compared with MC-SD data (dotted 
curves) for the L = 128 system. The peaks in the MC-SD data 
fall at slightly lower frequencies than the memory predictions. 
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FIG. 13. The T = 0.5J dynamical structure function 
obtained from the memory function calculations (solid curves, 
using L = 1024-2048) compared with MC-SD data (dotted 
curves) for the L = 128 system. 
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FIG. 15. Dynamic structure function for an odd integer 
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FIG. 17. Plot of difference process weights versus spin- 
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FIG. 19. Sum (+ signs) and difference (solid cir- 
cles) process weights versus frequencies Q± — ± tA> q _k 
for the L = 128 system, for an odd integer wavevector 
q = (5, 0)(27r/L), at T — 0.1J. Compare the locations of 
the features in Fig. [l(| Numbers in parenthesis are wavevec- 
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FIG. 20. Dynamic structure function at wavevectors 
q = (n/8, 7r/8) and q = (7r/4, 7t/4), from the memory calcula- 
tions (solid curves) compared to Monte Carlo Spin-Dynamics 
(data points) on L — 128 systems, for a range of higher tem- 
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